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Abstract 

We consider a fluid of d-dimensional spherical particles interacting via a pair potential <f>(r) which takes a finite value e if the 
two spheres are overlapped (r < a) and otherwise. This penetrable-sphere model has been proposed to describe the effective 
interaction of micelles in a solvent. We derive the structural and thermodynamic functions in the limit where the reduced 
temperature ksT/e and density pa d tend to infinity, their ratio being kept finite. The fluid exhibits a spinodal instability at 
a certain maximum scaled density where the correlation length diverges and a crystalline phase appears, even in the one- 
dimensional model. By using a simple free- volume theory for the solid phase of the model, the fluid-solid phase transition is 
located. 
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Most of the theoretical studies and numerical applica- 
tions of the theory of liquids in equilibrium is devoted to 
particles which interact according to unbounded spheri- 
cally symmetric pair potentials [1,2]. Atomic and molec- 
ular fluids have been usually modeled in this way and a 
vast effort was done during the second half of the past 
century in order to understand systems such as hard 
spheres, the square-well model, or the Lennard-Jones 
liquid. The development of the integral equations of the 
theory of liquids, the Percus-Yevick (PY) approxima- 
tion being among the most widely used in this context, 
and their solution for some simple models [1,2,3] were 
landmarks in the history of the theory of liquids. 

In the last decade, the properties of fluids of parti- 
cles interacting via bounded pair potentials have been 
the subject of an increasing interest, the Gaussian 
core model [4,5,6] and the penetrable-sphere model 
[5,6,7,8,9,10,11,12,13] being among the most popular- 
ones. The motivation for the study of fluids based upon 
this new class of interactions is two-fold. First, from a 
fundamental point of view, they are useful to unveil the 
weaknesses of standard integral equation theories and 
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other approximations whose validity has only been as- 
sessed from applications to unbounded potentials. More 
consistent closure approximations arise from these stud- 
ies [10,11]. From a more practical point of view, these 
models have also been proposed in order to understand 
the peculiar behavior of some colloidal systems, such 
as micelles in a solvent or star copolymer suspensions. 
The particles in these colloids are constituted by a small 
core surrounded by several attached polymeric arms. As 
a consequence of their structure, two or more of these 
particles allow a considerable degree of overlapping with 
a small energy cost [6]. An ultrasoft logarithmically 
divergent potential for short distances has also been 
proposed to describe the effective interaction between 
star polymers in good solvents [14]. These are a few 
particular cases of systems defining what is commonly 
known as "soft matter" , which has become an active 
field of research with many potential physical, chemical 
and engineering applications. 

The penetrable-sphere (PS) interaction potential is de- 
fined as 



(r) 




(1) 



This model was suggested by Marquest and Witten [15] 
in the late eighties as a simple theoretical approach to 
the explanation of the experimentally observed crystal- 
lization of copolymer mesophases, where a simple cu- 
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bic solid phase coexists with the disordered suspension. 
By arguments based on the internal energy alone, these 
authors claimed that, under the assumption of single- 
site occupancy, the stability of the simple cubic crystal 
was assured in a given domain of the phase diagram. 
On the other hand, density-functional theory [8,9] pre- 
dicts a freezing transition to fee solid phases with mul- 
tiply occupied lattice sites. The existence of clusters of 
overlapped particles (or "clumps" ) in the PS crystal and 
glass was already pointed out by Klein et al. [7], who 
also performed Monte Carlo simulations on the system. 
In the fluid phase, the standard integral equation the- 
ories (e.g., PY and HNC) are not very reliable in de- 
scribing the structure of the PS fluid, especially inside 
the core [8]. Other more sophisticated closures [10], as 
well as Rosenfeld's fundamental-measure theory [11] , are 
able to predict the correlations functions with high pre- 
cision but, on the other hand, the signature for a spa- 
tially ordered phase [a divergent structure factor S(k) 
for a finite wavenumber k] is not found with these meth- 
ods. Recently, a mixture of colloids and non-interacting 
polymer coils has been studied [12], where the colloid- 
colloid interaction is assumed to be that of hard spheres 
and the colloid-polymer interaction is described by the 
PS model. The inhomogeneous structure of penetrable 
spheres in a spherical pore has also been investigated 
[13]. 

In order to shed further light on the properties of the PS 
system, in this Letter we focus on the high-temperature, 
high-density region of the phase diagram. As we will see 
later, this asymptotic region is mathematically equiva- 
lent to taking e — > and a — > oo in a scaled way, so 
that we recover rigorous results first derived by Gates 
[16] and by Grewe and Klein [17] in the more general 
framework of a class of Kac potentials. At the level of 
the thermodynamic properties, only zeroth-order results 
were given in Ref. [17], but here we provide the explicit 
first-order corrections. Moreover, we combine the exact 
free energy of the PS fluid with a free-volume estimate 
for the free energy of the PS crystal to obtain the freez- 
ing and melting points. 



x w e/ksT, it is more convenient to work with x rather 
than with fc^T/e as the control parameter. 

Let us consider now the exact virial expansion of the 
cavity function y(r) = g(r) exp[(f>(r) / ksT], where g(r) is 
the radial distribution function [1]. In the diagrammatic 
representation of the virial expansion of y(r), each bond 
in a given diagram represents a Mayer function /pg (r) = 
xfns(r)- Therefore, 



y(r) 
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E 
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x m yW(r), (3) 



where the HS function yi"^ (r) is represented by the sum 
of diagrams having 2 root points (white circles) , n field 
points (black circles) and n + 1 + m bonds. Setting x — 
1 (zero-temperature limit), Eq. (3) becomes the virial 
series for hard spheres. On the other hand, in the high- 
temperature limit x — > with p — > oo and p — px finite, 
we get 



y(r) = l+xJ2p n y { n ) (r) + 0(x 2 ) = l + xw(r)+0(x 2 ), 

71=1 

(4) 

where the second equality defines the function w(r). The 
functions yn\r) are represented by chain diagrams, i.e., 



yM(\ ro -r n+1 \) = t±^- 



n n+l , . 

-•— o . (5) 



By application of the convolution theorem, the Fourier 
transform of the function w(r) is given by 



w(k) = p[Ms(k)} 2 l-Afas(fc) 



(6) 



where fns(k) = — (27rer) d//2 fc d l 2 Jd/2{ka) is the Fourier 
transform of the HS Mayer function, J u (z) being the 
Bessel function. Similarly, the structure factor S(k) and 
the direct correlation function c(r) are 



We start by observing that the Mayer function of the PS 
interaction is simply 



S(k)= 1-p/hsW +0(x) 



(7) 



/p S (r) = cxp [-4>(r)/k B T] - 1 = x/ HS (r), (2) 

where x = 1 — e~ e / fcsT is a parameter measuring the 
temperature of the system, /hsM = — 0(<r — r) is the 
Mayer function of a hard-sphere (HS) system (6 be- 
ing the Heaviside step function), ks is the Boltzmann 
constant and T is the absolute temperature. Obviously, 
the PS model includes the HS fluid (T -> or x -> 1) 
and the point-particle fluid (T — > oo or x — > 0) as spe- 
cial limits. The latter limit is assumed at finite num- 
ber density p. However, a non-trivial limit is obtained if 
x — ► and p — > oo keeping the product (scaled density) 
p = px finite. Although for high temperatures one has 



c(r)=xfn S {r)+0(x 2 ). (8) 

From Eqs. (4) and (8) it follows that the PY and HNC 
closures become exact to first order in x. The same hap- 
pens with any sensible approximation which retains the 
chain diagrams of the virial expansion. It is worth noting 
that the non-linear Debye-Huckel for the radial distri- 
bution function of a system of charged particles can be 
derived by neglecting all but the chain diagrams, which 
are most weakly connected (and hence the most strongly 
divergent) ones [1] . More in general, the chain diagrams 
determine the asymptotic long-distance behavior of the 
correlation functions [2] . 



2 



Table 1 

Contributions of zeroth- and first-order in x to th 
thermodynamic quantities. 



Quantity 


O(x ) 


0(x) 


a cx /k B T 


2 d ~ 1 rj 


2 d - 1 7 (a) 


p/pk B T 


l + 2 d - 1 rj 


2 d - 1 rjw(a) 


Scx/fcs 





2 d -^[w(a)-\] 




2 d rf 


2 d " 1 [t}-2- d w(0)] 




2 d - 1 rj 


2 d -' 2 [n-2 1 - d w(0)] 







±rjdw(0)/drj 


{dp/dp) T /k B T 


f + 2 d rj 


2 d - 1 d[rj 2 w{a)\/drj 



Besides the standard structure functions of the th 
liquids listed above we have found it useful in th 
to introduce an auxiliary function j(r) as the ( 
integral of w(r), namely w(r) = d^{r)/drj, whe 
pvd<J d is the scaled packing fraction, va = (7r/4) d / 2 
d/2) being the volume of a d-dimensional sphere 
diameter. In Fourier space, 



7(fc) = -rjfns(k) - v d a d In 1 - pfns(k) 



We have verified that the values of 7(7") at r = 
r = a satisfy the linear relation 



+ —7(0) 



1 

2' 



which is crucial to prove the thermodynamic cons: 
between the virial and energy routes to the equa 
state. By standard application of statistical-mecl 
formulas relating the correlation functions to th 
modynamic quantities [1], we have derived expr 
for the latter to first order in x in terms of the ve 
w(r) and 7 (r) at r = and r = a. The (excess) f 
ergy per particle a cx , the compressibility factor^/ 
the (excess) entropy per particle s ex , the (excess) 
ical potential /i cx , the (excess) internal energy p 
tide u ox , the (excess) specific heat c cx and the 
isothermal compressibility (dp/dp^/ksT are li; 
Table 1. The second and third columns give the : 
and first-order contributions, respectively, in th< 
expansions of those quantities in powers of x at cc 
p = px. The zeroth-order terms are linear functic 
(or rj), i.e., the virial expansion truncated after t 
ond virial coefficient becomes exact in the limit 
However, the first-order terms are highly nonlinea 
tions of the scaled density and so all the virial coef 
contribute. 

Gates [16] and Grewe and Klein [17] consid 
class of Kac potentials of the form <p(r) — e<j>' 
where <p*{r*~) is a non-negative bounded anc 
grable function and e oc a . They proved rigc 
that in the van der Waals limit a — > 00, the 
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Fig. 1. Plot of s e x/(e/T) (lower curves) and Ce X /(e/T) (upper 
curves) versus 77/770 for d = 1 (dotted line), d = 2 (dashed 
line) and d = 3 (solid line). 

correlation function adopts the mean-field expression 
c(r) = —<j>{r)/kBT and hence the structure factor is 



S(k) = [l + (p/k B T)<P(k) 

to the Kac class with (f>*(r*) 
van der Waals limit one has x 



. The PS potential belongs 

= 0(1 — r*) and so in the 
— > with rj fixed. It is then 



obvious that the van der Waals limit (e cx cr , a — > 00 
at fixed T and p) is equivalent to the high-temperature, 
high-density limit (at fixed e and cr) considered in this 
paper. In fact, Eqs. (7) and (8) are recovered from the 
more general results of Refs. [16] and [17], although the 
route followed here differs from theirs. Equations (7) 
and (8) were also found by Likos et al. [5] as a mean- 
field approximation in the limit p — > 00. However, as is 
apparent from our derivation of Eqs. (7) and (8), the 
key ingredient is the high-temperature assumption. The 
additional high-density assumption is only needed to 
depart from the trivial results corresponding to a gas 
of non-interacting particles. Comparison with Monte 
Carlo simulations shows that the mean-field structural 
functions behave very well even for x < 0.2 [5]. 

The results for the pressure, the internal energy and 
the isothermal compressibility to zeroth-order in x were 
already given by Grewe and Klein [17]. However, Eq. 
(10) and the third column of Table 1 are, to the best of 
our knowledge, new results. The contributions of 0{x) 
are especially relevant in the case of the excess entropy 
and specific heat, since the corresponding O(a; )-terms 
vanish. Figure 1 shows the high-temperature limit of 
s cx /(e/T) and c ex /(e/T) versus 77/770 (where 770 will be 
defined below). The nonlinear density-dependence of 
both quantities is clearly observed. 

As analyzed by Grewe and Klein [17,18], the PS fluid 
presents a spinodal point (Kirkwood instability) at a cer- 
tain scaled density. A simple examination of Eqs. (6), 
(7) and (9) shows that these quantities are well defined 
for every real wavenumbcr k if and only if the maxi- 
mum value of fns{k) is smaller than l/p. In general, the 
absolute maximum / max of fas(k) occurs at k , where 



3 



Table 2 

Values of the HS close-packing fraction rj cp , the wavenumber 
feo, the nearest-neighbor distance ro, the (scaled) spinodal 
instability packing fraction 77b, the (scaled) freezing packing 
fraction rj{, the (scaled) packing fraction 5) ms at the condi- 
tion of marginal stability and the (scaled) melting packing 
fraction rjm. 



tion w(r) behaves as w(r) w (1 — 77/770) 1 ^ 2 to{r), where 
the scaling function uj(r) is given by 



Lu(r) 



(27T) 
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V2a 



(fc 0( r)^ 1 (fcor)- d/2+1 J (i/2 _ 1 (fcor)- 
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1.40 


2.30 


1.00 


1.00 


1.00 


2 


v^Tr/e ~ 


0.91 


5.14 


1.37 


1.89 
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0.95 


1.03 
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0.74 


5.76 


1.34 


1.45 


0.62 


0.69 


0.80 


4 


tt 2 /16 ~ 


0.62 


6.38 


1.32 


1.07 


0.36 


0.41 
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V / 2tt 2 /30 : 


- 0.47 
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0.22 


0.26 


0.33 
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Fig. 2. Plot of w(r) and S(k) (see inset) at 77/770 = 0.1 (dotted 
line), 0.5 (dashed line) and 0.9 (solid line) for d = 1. 

ko<7 is the first zero of J<j/2+i( z )- Therefore, there ex- 
ists an upper bound po = 1/ / max to the scaled density, 
such that the structure factor becomes divergent at the 
wavenumber k = ko when p — > p . The values of ko 
and 770 = Po v d& d for d = 1-5 are displayed in the third 
and fifth columns of Table 2, respectively. As an illus- 
tration, in Fig. 2 we have plotted the functions w(r) and 
S(k) corresponding to the one-dimensional PS fluid for 
77 = 0.1?7o, O.5770 and 0.9t7o. 

How do the structural and thermodynamic functions 
behave as the scaled packing fraction 77 approaches its 
upper bound 77b from below? Let us denote by k = 
±k(t7) ± 19(77) (with the convention n, q > 0) the four 
zeroes of 1 — jo/us (k) closest to the real axis in the com- 
plex fc-plane. These values are responsible for the asymp- 
totic behavior of w(r) [and hence of g(r)] for long dis- 
tances: k is the wavenumber of the oscillations, while q 
is the damping coefficient, i.e., £ = q^ 1 is the correla- 
tion length. A straightforward asymptotic analysis yields 
ko~ ks k cr — [(d+ l)/3fco cr ](l — v/Vo) and 



(12) 

So, the first-order contribution w(r) to the cavity and 
radial distribution functions diverges as 77 tends to the 
maximum value 770- On the other hand, the auxiliary 
function 7(7") remains finite in that limit. According to 
Table 1, the first-order coefficients in x of the pressure, 
the entropy, the chemical potential and the internal en- 
ergy diverge as (1 — rj/rjo) -1 / 2 when rj — > 770, while the 
coefficients of the specific heat and the isothermal com- 
pressibility diverge as (1 — rf/rfo)~ 3 ^ 2 - Equation (12) im- 
plies that uj(t) ~ (fcor) - ( d-1 )/ 2 cos[fe 7' - (d- for 
large kor, so the correlations at 77 = rjo oscillate with 
a wavenumber fco, the amplitude decaying algebraically. 
The first maximum of u(r) (apart from the one at r = 0) 
occurs at a value r = r such that kor is the second 
zero of Jdii{z)- The quantity ro represents the distance 



>d/2 

between nearest neighbors at 77 
2 for d = 1-5. 



770 • It is given in Table 



On physical grounds, it is expected that the freezing 
transition from the stable fluid to the stable solid occurs 
at a scaled density % smaller than the value 770 at the 
spinodal instability. We have estimated the values of the 
scaled packing fraction at freezing (%), at melting (77m) 
and at the point of marginal stability (^s), by using for 
the PS solid phase a simple free-volume theory based 
on the one for the HS solid [8,20]. In the basic picture 
of the PS solid [6,7,8], the lattice sites are occupied by 
"clusters" of spheres that overlap each other. A parti- 
cle inside one of these clusters performs random motions 
but excursions beyond the lattice distance are forbidden 
because of the large energy cost associated with overlap- 
pings with particles in the neighbor clusters. In the high- 
temperature limit these clusters are expected to contain 
typically a large number of particles. We will suppose 
that, on average, this number scales as a/x (where a 
is a density-dependent parameter to be determined), so 
rj/ia/x) = rj/a is the packing fraction of the clusters. 
Consequently, the excess internal energy per particle of 
the PS solid is w|° lld /ksT — a/2. Under the assumption 
that every cluster behaves as a hard-core particle with 
a free volume 7jf roo (77/a), where 



V free (v) = 



e/a ^(qa)- 1 ^ 2~ 1 / 2 {1 -77/770)" 



1/2 



(11) 



2 d 
P 



1 - 



Vcp 



1/d' 



(13) 



which implies that the correlation length diverges with 
a critical exponent v = \ as the density approaches the 
maximum density [19]. Moreover, the correlation func- 



is an estimate of the free volume of a hard sphere in 
a crystal with packing fraction 77 [8,20] and 77 cp is the 
HS close-packing fraction, we have estimated the excess 
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Fig. 3. Excess free energy per particle in the 
three-dimensional PS solid, al^ ld /kBT (solid line), and PS 
fluid, aex ld /fesT (dashed line), in the high-temperature 
limit. The excess internal energy, u s ^ lA /UbT (dotted line), 
and the excess entropy, s| x h As (dashed-dotted line), of 
the PS solid are also plotted. The shaded area represents 
the fluid-solid coexistence region. 

entropy per particle of the PS solid as 



olid 



k B T 



= -din 



1 - 



1/dC 



-din 2. 



(14) 



Given a value of 77, the parameter a is determined by 



minimizing the free energy a; 



solid 

cx 



.solid 



+ Ts- 



olid 



111 



this way we find that a(fj) is the solution of the algebraic 
equation rj/rj cp = a d+1 / (2 + a) d , so that 



knT 2 



din 



a(rf) 



(15) 



The pressure and the excess chemical potential of the 
high-temperature PS solid are then p solld j pk B T = 
1 + a(rf)/2, v s c f d (v)/k B T = a(rj) + din [1/2 + a(r?)/4]. 
In Fig. 3 we have plotted the excess free energy per 
particle of the solid and the fluid phases for d = 3. 
The excess internal energy and the excess entropy 
of the solid are also plotted. For the fluid we have 
used the zeroth-order approximation given in Table 
1, so a^ id {fj)/k B T » u^ d (rj)/k B T « 2 rf - 1 f? and 
s cx ld (v)/kB ~ 0. While the solid has a smaller entropy 
than the fluid, it requires less internal energy. As the 
density increases, the latter effect dominates over the 
former and the solid becomes more stable than the 
fluid. The (scaled) density of marginal stability f} ms is 
determined by the condition a^° hd = a^ ld . The (scaled) 
freezing and melting densities are obtained from the 
equality of the pressure and the chemical potential in 



bothfases:p Huld (7? f Hp solid (5 m ),^ x uid (??f) = /4x Ud (? m )- 
The values of fj m , rj ms and rjf are given in Table 2 for 
d = 1-5. We observe that our heuristic estimates for the 
characteristic densities of the fluid-solid transition are 
typically less than half the upper bound density 770- 

The one-dimensional (ID) case deserves some special 
comments. For d = 1, the curves representing the free 



energies of the solid and the fluid do not cross, but "kiss" 
each other at rj = 1 (i.e., they have the same value and 
slope at rj = 1). In fact, a|° lld (?7) < a^ ld (rj) not only for 
rj > 1 but also for rj < 1. This "exaggerated" stability of 
the ID PS solid for small densities is obviouly an artifact 
of the heuristic free-volume theory we have employed. 
Nevertheless, our free-volume theory contains the basic 
ingredients explaining that the high-temperature ID PS 
solid becomes more stable than the fluid at sufficiently 
large densities rj. This is further confirmed by the rig- 
orous existence of the spinodal instability at 770 — 2.30, 
where the metastable fluid ceases to exist and a continu- 
ous transition to a crystalline solid takes place. We must 
emphasize that the fluid-solid phase transition in the ID 
PS model is not forbidden by van Hove's theorem [21] 
because one of the hypothesis of the theorem, requiring 
the interaction potential to include an infinite repulsive 
core that does not allow full overlaps between particles, 
is not fulfilled [22]. Thus, the PS model provides one of 
the rare examples of one-dimensional models exhibiting 
phase transitions [22] . Whether the fluid-solid transition 
occurs near rj = 1 or whether it does not present a den- 
sity jump, as indicated in Table 2, are questions that 
need further theoretical and simulational work before 
being satisfactorily elucidated. 



In this paper we have focused on the high-temperature 
domain (x — ► 0) of the d-dimensional PS model. When 
the exact diagrammatic expansion of the cavity func- 
tion is considered, it turns out that x acts as an order- 
ing parameter, so that the first-order term contains only 
chain diagrams, which can be easily resummed, yielding 
mean-field expressions. If the unsealed packing fraction 
rj is kept finite, one arrives at the trivial case of non- 
interacting particles. If, on the other hand, one explores 
the high-density regime rj = rj/x, much more interesting 
results appear. To zeroth-order in x [17] the thermody- 
namic quantities are described by the second virial coef- 
ficient alone, but the first-order corrections exhibit a rich 
nonlinear dependence on rj. The fluid presents a spinodal 
instability at the upper bound density rjo(x) = rjo/x, 
but this is preempted by a first-order phase transition 
to the solid at the freezing density rj[(x) — rjijx. It 
seems natural to expect that a similar situation ap- 
plies for finite and low temperatures, except that the IE- 
dependence of 770(2;) and rjiix) is more complicated than 
in the high-temperature limit. From that point of view, 
it can be conjectured that the zero-temperature limit 
(where the PS model reduces to the HS model) of the up- 
per bound density is Xm\ x ^,\rj^{x) = 77 cp . Analogously, 
linij;^! rjf(x) = rj^ s , where ?y^ s is the freezing packing 
fraction of the HS fluid. Since rj cp < rj (cf. Table 2) and 
rjf s < fjf, it seems plausible that the products xrj (x) and 
xrj{(x) are smoothly decreasing functions of x bounded 
between the values 770 and rjf, respectively, at x = and 
the values ri cp and T]f s , respectively, at x = 1. This sug- 
gests the possibility of constructing a simple approxi- 
mate theory for the PS model spanning the whole tem- 
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perature spectrum, by interpolating between the suc- 
cessful PY theory for hard spheres (zero temperature) 
and the mean- field-theory results (high temperatures). 
We are currently working along these lines and further 
results will be published elsewhere. 
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